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Abstract 



Our contribution concerns with the numerical solution of the 3D general 
relativistic hydrodynamical system of equations within the framework of the 
■ {3 + 1} formalism. We summarize the theoretical ingredients which are 

■^j- ' necessary in order to build up a numerical scheme based on the solution 

£C) , of local Riemann problems. Hence, the full spectral decomposition of the 

1 Jacobian matrices of the system, i.e., the eigenvalues and the right and 

left eigenvectors, is explicitly shown. An alternative approach consists in 
using any of the special relativistic Riemann solvers recently developed for 
describing the evolution of special relativistic flows. Our proposal relies 
on a local change of coordinates in terms of which the spacetime metric 
is locally Minkowskian and permits an accurate description of numerical 
Oh general relativistic hydrodynamics. 
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1 Introduction 

Astrophysical scenarios involving relativistic flows have drawn the attention and 



efforts of many researchers since the pioneering studies of May and White 1 16 and 



Wilson |23|]. Relativistic jets, accretion onto compact objects (in X-ray binaries 



CO 
> 

• i-H 

X. 

or in the inner regions of active galactic nuclei), stellar core collapse, coalescing 
compact binaries (neutron star and/or black holes) and recent models of formation 
of Gamma-ray bursts (GRBs) are examples of systems in which the evolution of 
matter is described within the frame of the theory of relativity (special or general) . 

Since 1991 |l4| the use of Riemann solvers in relativistic hydrodynamics has 
proved successful in handling complex flows, with high Lorentz factors and strong 
shocks, superseding traditional methods which failed to describe ultrarelativistic 
flows |l7| . Exploiting the hyperbolic and conservative character of the relativistic 
hydrodynamical equations, we proposed how to extend modern high-resolution 
shock- capturing (HRSC) methods to the relativistic case, first in one-dimensional 
calculations |14j, and, later on, we extended them to multidimensional special 
relativistic [j7| and multidimensional general relativistic hydrodynamics Qj. We 
made use of a linearized Riemann solver based on the spectral decomposition of 
the Jacobian matrices of the system. 
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Unlike the case of classical fluid dynamics the use of HRSC techniques in the 
frame of relativistic fluid dynamics is very recent and has yet to cover the full set of 
possible applications. Up to now, the most interesting astrophysical applications 
have involved the simulation of extragalactic relativistic jets (see [El and for, 
respectively, relativistic 3D-hydro and 2D-magnetohydro calculations). Recently, 
some studies on the morphology of accreting flows onto moving black holes have 
being carried out (see |8| and references therein) using a multidimensional general 
relativistic hydrocode. A very promising application of HRSC techniques in the 
frame of general relativistic magnetohydrodynamics has been used recently to 
simulate the formation of relativistic jets from black holes magnetized accretion 
disks @. 

At present, to develop robust and accurate general relativistic hydrocodes is a 
challenge in the field of Relativistic Astrophysics. A general relativistic hydrocode 
is a useful research tool for studying flows which evolve in a background spacctimc. 
Furthermore, when appropriately coupled with Einstein equations, such a general 
relativistic hydrocode is crucial to model the evolution of matter in a dynamical 
spacetime. The coupling between geometry and matter arises through the sources 
of the corresponding system of equations. Such a marriage between numerical 
relativity and numerical relativistic hydrodynamics could be useful, for example, 
to analyze the dynamics (and the physics) of coalescing compact binaries. These 
are one of the most promising sources of gravitational radiation to be detected 
by the near future Earth-based laser-interferometer observatories of gravitational 
waves. 

2 The equations of general relativistic hydrody- 
namics as a hyperbolic system of conservation 
laws 

The evolution of a relativistic fluid is governed by a system of equations which 
summarize local conservation laws: the local conservation of baryon number, V • 
J = 0, and the local conservation of energy-momentum, V ■ T = (V- stands for 
the covariant divergence). 

If {dt,di} define the coordinate basis of 4- vectors which are tangents to the 
corresponding coordinate curves, then, the current of rest-mass, J, and the energy- 
momentum tensor, T, for a perfect fluid, have the components: J M = pu^ , and 
T^ v = phu^u v +pgV, respectively, p being the rest-mass density, p the pressure 
and h the specific enthalpy, defined by h = I + s + p/ p, where e is the specific 
internal energy, is the four- velocity of the fluid and g^ v defines the metric of the 
spacetime M. where the fluid evolves. As usually, Greek (Latin) indices run from 
to 3 (1 to 3) - or, alternatively, they stand for the general coordinates {t, x, y, z} 
({x, y, z}) - and the system of units is the so-called geometrized (c = G = 1). 

An equation of state p = p(p, e) closes, as usual, the system. Accordingly, 
the local sound velocity c s satisfies: hc\ = x + (p/p 2 )n, with \ — dp/dp\ e and 
k = dp/de\ p . 
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Following jl , let M. be a general spacetime, described by the four dimensional 
metric tensor g^ v . According to the {3 + 1} formalism, the metric is split into the 
objects a (lapse), f3 l (shift) and 7^, keeping the line element in the form: 

ds 2 = -(a 2 - fiip^dt 2 + 2f3 t dx'dt + ltj dx l dx j (1) 

If n is a unit timelike vector field normal to the spacelike hypersurfaces T, t (t 
= const.), then, by definition of a and (i l is: d t = an + (3 l di, with n • di = 0, 
Vi. Observers, Oe, at rest in the slice £<, i.e., those having n as four-velocity 
(Eulerian observers) , measure the following velocity of the fluid 

cm 1 a 

where W = — (u • n) = an*, the Lorentz factor, satisfies W = (1 — v 2 ) -1 / 2 with 

V 2 = V^ 1 (Vi = 7yV). 

Let us define a basis adapted to the observer Oe, e^) = {n, di}, and the 
following five four-vector fields {J, T n, T ■ d\, T • <9 2 , T • 93}. Hence, the above 
system of equations of general relativistic hydrodynamics (GRH) can be written 

V • A = s, (3) 

where A denotes any of the above 5 vector fields, and s is the corresponding source 
term. 

The set of conserved variables gathers those quantities which are directly mea- 
sured by Oe, i.e., the rest-mass density (D), the momentum density in the j- 
direction (Sj ) and the total energy density (E) . In terms of the ■primitive variables 
w = (p, Vi, e) (vi = ji-jvi) they are 

D = pW , Sj = phW 2 v 3 , E = phW 2 - p (4) 

Taking all the above relations together, the fundamental system to be consid- 
ered for numerical applications is 



1 { dy/jF° (yr) dy/=gF*(w) 



s(w) (5) 



\J—g \ dx° dx l 
where the quantities F Q (w) are 

F°(w) = (A5 i) r) (6) 

F'(w) = (D (v* - £) , Sj (v* - £) + p*J,r („« - t) + V) (7) 
and the corresponding sources s(w) are 

s(w) = (0, T- - 1*^) , a (t*^ - T-I^) ) (8) 

r being r = E — D, and g = det(g M „) is such that y/—g = a^fj (7 = det(7y)) 
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2.1 Linearized Riemann Solvers and Characteristic fields 



Modern HRSC schemes use the characteristic structure of the hyperbolic system 
of conservation laws. In many Godunov-type schemes, the characteristic structure 
is used to compute either an exact or an approximate solution to a sequence of 
Riemann problems at each cell interface. In characteristic based methods the 
characteristic structure is used to compute the local characteristic fields, which 
define the directions along which the characteristic variables propagate. In both 
these approaches, the characteristic decomposition of the Jacobian matrices of the 
nonlinear system of conservation laws is important, not only because it is one 
of the key ingredients in the design of the numerical flux at the interfaces, but 
because experience has shown that it facilitates a robust upgrading of the order 
of a numerical scheme. 

The three 5 x 5- Jacobian matrices B % associated to system (||) are 



B' 



a 



The eigenvalues of B x are 



A = av x - f3 x (triple) 



(9) 



(10) 



which defines the material waves, and two other, A±, associated with the acoustic 
waves 



A± = 



o 



1 — v 2 c'i 



K(i-c2)± 



±c s y/(l - w 2 )[7^(l - v 2 c 2 s ) - v x v x (l - C 2)]| - f3 x 
A complete set of right- eigenvectors is 



(11) 



hW 



1 



v x - A| 

^yXX — . y x 



hWv„ 



hWv 7 



hW(pi xx — v x v x ) 



v x Al 



ro,i = 



K_ 

hW 



hW 
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1"0, 2 



ro,3 



WVy 

h ( lxy + 2W 2 v x v y ) 

h ( ly y + 2W 2 VyVy) 
h ( lz y + 2W 2 V Z Vy) 

Wv y (2hW-l) 
The corresponding set of left- eigenvectors is 

1 M = (h - W, Wv x , Wv y , Wv z , -W) 

K, — 1 



Wv z 

h ( lxz + 2W 2 v x v z ) 
h ( lyz + 2W 2 v y v z ) 
h { lzz + 2W 2 v z v z ) 
Wv z (2hW-l) 



l 0,2 



-JzzVy 
V X {lzzVy -ly Z V z ) 

7^(1 - v x v x ) + i xz v z v x 

-7^(1 - - IxzVyV* 

-IzzVy + lyzVz 



1(3,3 — 



lyy v z + JzyVy 

V X {lyyV z - lzy Vy) 

-7^(1 - - IxyVzV 3 - 

-7i/i/«« + 7*i/«i/ 



h 2 



r xx (l - /C4^) + (2/C - 1)V£(WV£ - T xx v x ) 
T xy {\ - KA X ± ) + (2/C - 1)V£(WV£ - r^v*) 

r„(i - ACi|) + (2/c - \)v x ± {w 2 v^ - r,/) 

L {l-K,)[-~(v x + Vl{W 2 £-T xx )]-KW 2 VU 
where the following auxiliary quantities and relations have been introduced: 

A l ± = A±+/3 4 , A = A/a , p l = f3 l /a 



AC = 



, K — K,/p 



(12) 
(13) 
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six \ r 

C± = v x - V. 




ryxx — y x 



(14) 




(15) 



A = h 3 



(16) 



i = T xx - jv x v 



.X 



j Txx — lyy^fzz 7;, 



(17) 



Symmetry arguments allow one to obtain the spectral decomposition in the 
other spatial directions. The corresponding expressions in Special Relativity j^] 
are easily covered. The above full spectral decomposition provides the user with 
the technical ingredients needed to develop state-of-the-art, upwind-based HRSC 
codes for numerical relativistic hydrodynamics. In 3D general-relativistic applica- 
tions, and depending on the particular Riemann solver or flux formula used, the 
knowledge of the analytical values of the left-eigenvectors could be crucial in the 
efficiency of the code (see [||). 

3 Special Relativistic Riemann Solvers in General 
Relativistic Hydrodynamics 

Up to now, only a small number of papers have considered the extension of HRSC 
methods to GRH using linearized Riemann Solvers or flux formulae (|lj] and j2^] 
for ID problems, || and G3|). or deriving explicitly an extension of Roe's 
Riemann solver to GRH (||). 

In |H| we have answered the following basic question (suggested in || and jl3| ) : 
Is it possible to obtain a general procedure that allows one to take advantage 
of Special Relativistic Riemann Solvers (SRRS) to generate numerical solutions 
describing the evolution of relativistic flows in strong gravitational fields? 

The affirmative answer relies on a local change of coordinates, at each numerical 
interface, in terms of which the spacetime metric is locally Minkowskian. Our 
procedure, hence, follows analogous trends to those used in classical fluid dynamics 
to solve Riemann problems in general curvilinear coordinates. 

Let us consider the integral form of the system of equations (|^) on a four- 
dimensional volume f2, with three-dimensional boundary dfl, and apply Gauss 
theorem to obtain the corresponding balance equation 



For numerical applications, we choose volume f2 as the one bounded by the 
coordinate hypersurfaces {S^q, T, x a + ^ x a}. Hence, the time variation of the mean 
value of A , 




(18) 




(19) 
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within the spatial volume 



fX 1 +Ax 1 ,-x 2 +Ax 2 px 3 +Ax 3 

AV = / / / ^gdx 1 dx 2 dx 3 , (20) 



can be obtained from 



(A°AV) t +At = (A°AV) t + [ sdn 



A • d S + / A-d J S 

il+Ax 1 



A-d 3 S+ / A-d 3 S 



/ A ■ d S + / A ■ d X I . (21) 

In order to advance in time, the volume and surface integrals on the right hand 
side have to be evaluated. We have applied HRSC to calculate the A vector fields 
by solving local Ricmann problems combined with monotonized cell reconstruction 
techniques. 

According to the Equivalence Principle, physical laws in a local inertial frame 
of a curved spacetime have the same form as in special relativity (see, e.g., Schutz 
[p2j). Locally flat (or geodesic) systems of coordinates, in which the metric is 
brought into the Minkowskian form up to second order terms, are the realization 
of these local inertial frames. However, whereas the coordinate transformation 
leading to locally flat coordinates involves second order terms, locally Minkows- 
kian coordinates are obtained by a linear transformation. This fact is of crucial 
importance when exploiting the self-similar character of the solution of the Rie- 
mann problems set up across the coordinate surfaces. 

Hence, we propose to perform a coordinate transformation to locally Minkows- 
kian coordinates at each numerical interface assuming that the solution of the 
Riemann problem is one in special relativity and planar symmetry. This last as- 
sumption is equivalent to the approach followed in classical fluid dynamics, when 
using the solution of Riemann problems in slab symmetry for problems in cylin- 
drical or spherical coordinates, which breaks down near the singular points (e.g. 
the polar axis in cylindrical coordinates). Analogously to classical fluid dynam- 
ics, the numerical error will depend on the magnitude of the Christoffel symbols, 
which might be large whenever huge gradients or large temporal variations of the 
gravitational field are present. Finer grids and improved time advancing methods 
will be required in those regions. 

In the rest of this section we will focus on the evaluation of the first flux integral 



in Eq. (21). 



d 3 S = / A 1 y /=g~dx°dx 2 dx 3 (22) 
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To begin, we will express the integral on a basis with eg = n and ej forming an 
orthonormal basis in the plane orthogonal to n with ej normal to the surface T, x i 
and and eg tangent to that surface. The vectors of this basis verify e<j -e^ = rj & ^ 
with r\ & o the Minkowski metric (in the following, caret subscripts will refer to 
vector components in this basis). 

Denoting by Xq the coordinates at the center of the interface at time i, we 
introduce the following locally Minkowskian coordinate system 

x & = M & Jx a -x«), (23) 

where the matrix M£ is given by d a = M"ea, calculated at Xq. In this system 
of coordinates the equations of general relativistic hydrodynamics transform into 
the equations of special relativistic hydrodynamics, in Cartesian coordinates, but 
with non-zero sources, and the flux integral (|2^) reads 

/ (A 1 - —A°)\/^g~dx°dx^dx^ (24) 

with y/— g = 1 + 0(x a ), where we have taken into account that, in the coordinates 

x a , Sj.1 is described by the equation x^ — ^x° — (with j3 l — Mlf3 l ), where the 
metric elements (3 1 and a are calculated at Xq. Therefore, this surface is not at 
rest but moves with speed (3 1 /a. 

At this point, all the theoretical work on SRRS developed in recent years, can 
be exploited. The quantity in parenthesis in ( |2~i| ) represents the numerical flux 
across S^i , which can, now, be calculated by solving the special relativistic Rie- 
mann problem defined with the values at the two sides of T, x i of two independent 
thermodynamical variables (namely, the rest mass density p and the specific inter- 
nal energy e) and the components of the velocity in the orthonormal spatial basis 
v % (v l — M^v 1 ). Although most linearized Ricmann solvers provide the numerical 
fluxes for surfaces at rest, it is easy to apply them to moving surfaces, relying on 
the conservative and hyperbolic character of the system of equations (as in jn| ) . 

Once the Riemann problem has been solved, by means of any linearized or 
exact SRRS, we can take advantage of the self-similar character of the solution of 
the Riemann problem, which makes it constant on the surface T, x i simplifying the 
calculation of the above integral enormously ( pi| ) : 

/ A • d 3 T, = (A^ — —A°)* [ y^dx^dx^dx* (25) 

where the superscript (*) stands for the value on S a i obtained from the solution 
of the Riemann problem. The integral in the right hand side of (Q) is the area of 
the surface T, x i and can be expressed in terms of the original coordinates as 

/ xfi^./^gdx^dx^dx 5 (26) 
which can be evaluated for a given metric. 
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Finally, notice that the numerical fluxes defined in (|2J) correspond to the vector 
fields A = {J, T • n, T • ej, T • eg, T eg}. Thus the additional relation 

T-d i = M~i(T-e- j ) (27) 

has to be used for the momentum equations. 

The interested reader can address reference [^0) for details on the testing and 
calibration of our procedure. The additional computational cost of the approach is 
completely negligible. The procedure has a large potentiality and can be applied 
to other systems of conservation laws, such as magneto- hydrodynamics (MHD), 
making possible to solve the general relativistic MHD equations using the corre- 
sponding Ricmann solvers developed for the special relativistic case. 



4 Conclusions 

An appropriate conservative formulation for the equations, together with the 
knowledge of the characteristic fields associated to the system, define the starting 
point for using HRSC schemes. The spectral decomposition of the Jacobian matri- 
ces, corresponding to the fluxes in each spatial direction, is used in the numerical 
flux computation and, moreover, it is potentially interesting in allowing an exten- 
sive range of application of HRSC methods with different approximate Riemann 
solvers or flux formulae. 

The procedure outlined in Section §3 is -from the computational point of view 
- very cheap, since it involves a linear change of coordinates. It has a large 
potentiality and can be applied to other systems of conservation laws, as magneto- 
hydrodynamics, giving a very useful numerical tool to solve the general relativistic 
MHD equations using the corresponding Riemann solvers developed for the special 
relativistic case. In particular, it is possible to use the exact solution of the special 
relativistic Riemann problem jL5l, flljJ. 

The astrophysical applications foreseen in the present and near future include 
the study of jet formation, multidimensional stellar core collapse, gamma-ray 
bursts and the coalescence of compact binaries. HRSC methods can, without 
doubt, be used successfully to tackle these scenarios, and acquire the prestige they 
have already earned in the simulation of relativistic jets and accretion flows around 
compact objects. 
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